% ex02            Example 2: variants of 9/7 wavelet filters
% which is the same as "79 filter"

%----------------------------------------------------------------------
% Copyright (c) 2010.  Karl Skretting.  All rights reserved.
% Stavanger University, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  15.03.2010  KS: m-file made
%----------------------------------------------------------------------

A = double(imread('lena.bmp'))-128;
[M,N] = size(A);
% define a function for Coding Gain
funCG = @(x) mean(x(:))/geomean(x(:));  % x should be variance of subbands
level = 3;  Ms = 2^level; Ns = Ms;

if exist('myim2col','file') == 2
    disp('   The ''j97'' implementation in myim2col / mylwt2 (mirroring at edges):');
    disp('JPEG-2000 9/7 wavelet and scaling which emphazis DC.');
    Y = myim2col(A, 't','j97', 's',[Ms,Ns]);
    sigma2 = var(Y,0,2);  % variance for each row (subband)
    disp(['Estimated coding gain for ',int2str(Ms),'x',int2str(Ns),...
        ' ''j97'' wavelet is ',num2str(funCG(sigma2))]);
    disp(['Range of Y is ',num2str(min(Y(:))),' to ',num2str(max(Y(:)))]);
    Ar = mycol2im(Y, 't','j97', 'imsize',[M,N]);
    if (norm(A-Ar,'fro') < 1e-6)
        disp('restored image O K ');
    else
        disp('restored image N O T   O K');
    end
end

if exist('myim2col','file') == 2
    disp('   The ''m97'' implementation in myim2col / mylwt2 (mirroring at edges) :');
    disp('JPEG-2000 9/7 wavelet and scaling which equal energy in DC and HP.');
    Y = myim2col(A, 't','m97', 's',[Ms,Ns]);
    sigma2 = var(Y,0,2);  % variance for each row (subband)
    disp(['Estimated coding gain for ',int2str(Ms),'x',int2str(Ns),...
        ' ''m97'' wavelet is ',num2str(funCG(sigma2))]);
    disp(['Range of Y is ',num2str(min(Y(:))),' to ',num2str(max(Y(:)))]);
    Ar = mycol2im(Y, 't','m97', 'imsize',[M,N]);
    if (norm(A-Ar,'fro') < 1e-6)
        disp('restored image O K ');
    else
        disp('restored image N O T   O K');
    end
end

if (exist('getwave','file') == 2) && (exist('myim2col','file') == 2)
    disp('   The ''db79a'' as a matrix from getwave.m used in myim2col :');
    G = getwave('db79a',level);
    Y = myim2col(A, 't',G, 's',[Ms,Ns]);
    sigma2 = var(Y,0,2);  % variance for each row (subband)
    disp(['Estimated coding gain for  getwave( ''db79a'', level ) ',...
        ' transform is ',num2str(funCG(sigma2))]);
    disp(['Range of Y is ',num2str(min(Y(:))),' to ',num2str(max(Y(:)))]);
    disp('Synthesis uses  GetWave( ''db79s'', level )  ');
    F = getwave('db79s',level);
    Ar = mycol2im(Y, 't',F, 'imsize',[M,N]);
    if (norm(A-Ar,'fro') < 1e-6)
        disp('restored image O K ');
    else
        disp('restored image N O T   O K');
    end
end

if (exist('Decom2D','file') == 2) 
    % Decom2D(2, ...) uses getwave and Ttimes. Result is exactly as above
    disp('   Analysis is:  Decom2D(2, A, ''db79'', level)  which uses getwave');
    Y = Decom2D(2, A, 'db79', level);
    Y = im2col(Y, [Ms,Ns], 'distinct');
    sigma2 = var(Y,0,2);  % variance for each row (subband)
    disp(['Estimated coding gain is ',num2str(funCG(sigma2))]);
    disp(['Range of Y is ',num2str(min(Y(:))),' to ',num2str(max(Y(:)))]);
    disp('Synthesis is:  Decom2D(-2, Y, ''db79'', level)');
    Y = col2im(Y, [Ms,Ns], [M,N], 'distinct');
    % Y = mycol2im(Y, 's',[Ms,Ns], 'imsize',[M,N]); % OR myim2col with no transform!
    Ar = Decom2D(-2, Y, 'db79', level);
    if (norm(A-Ar,'fro') < 1e-6)
        disp('restored image O K ');
    else
        disp('restored image N O T   O K');
    end
end

if (exist('lwt2','file') == 2) && (exist('im2col','file') == 2)
    % this is like Decom2D with transform 4
    disp('   Matlab lwt2 lifting implementation (with zero padding)');
    disp('JPEG-2000 9/7 wavelet and scaling which emphazis DC.');
    LS = {'d', -1.586134342059924*[1,1], 1; ...
        'p', -0.052980118572961*[1,1], 0; ...
        'd',  0.882911075530934*[1,1], 1; ...
        'p',  0.443506852043971*[1,1], 0; ...
        1.230174104914001, 1/1.230174104914001, []};  % as in JPEG-2000
    %   1.1496, 1/1.1496, []};  % more equal distribution of energy
    Y = lwt2(A, LS, level);
    Y = im2col(Y, [Ms,Ns], 'distinct');
    sigma2 = var(Y,0,2);  % variance for each row (subband)
    disp(['Estimated coding gain is ',num2str(funCG(sigma2))]);
    disp(['Range of Y is ',num2str(min(Y(:))),' to ',num2str(max(Y(:)))]);
    Y = col2im(Y, [Ms,Ns], [M,N], 'distinct');
    Ar = ilwt2(Y, LS, level);
    if (norm(A-Ar,'fro') < 1e-6)
        disp('restored image O K ');
    else
        disp('restored image N O T   O K');
    end
    disp('Note that for lwt2 min and max are like ''m97'' (uses mylwt2 with mirroring),');
    disp('these two methods gives the same result, except on the edges.');
end
%
disp(' ');